The objectives of this notebook are:
library(AnnotationDbi)
library(org.Hs.eg.db)
library(Seurat)
library(Matrix)
library(tidyverse)
# Paths
path_to_obj <- here::here("scRNA-seq/results/R_objects/seurat_integrated_with_doublets.rds")
path_to_doubl_annot <- here::here("scRNA-seq/2-QC/5-batch_effect_correction/1-doublet_exclusion/tmp/doublet_final_annotations_and_pDNN.rds")
path_to_save <- here::here("scRNA-seq/results/R_objects/seurat_without_doublets_missing_reintegration.rds")
# Functions
source(here::here("scRNA-seq/bin/utils.R"))
tonsil <- readRDS(path_to_obj)
doublet_annot <- readRDS(path_to_doubl_annot)
if (all(rownames(doublet_annot) == colnames(tonsil))) {
tonsil$pDNN_hashing <- doublet_annot$pDNN_hashing
tonsil$pDNN_scrublet <- doublet_annot$pDNN_scrublet
tonsil$pDNN_union <- doublet_annot$pDNN_union
} else {
warning("barcodes are not equal")
}
pDNN_vars <- c("pDNN_hashing", "pDNN_scrublet", "pDNN_union")
pDNN_ggs <- purrr::map(pDNN_vars, function(x) {
p <- feature_plot_doublets(seurat_obj = tonsil, feature = x)
p
})
names(pDNN_ggs) <- pDNN_vars
pDNN_ggs
## $pDNN_hashing
##
## $pDNN_scrublet
##
## $pDNN_union
doublet_vars <- c("HTO_classification.global", "scrublet_predicted_doublet",
"has_high_lib_size")
doublet_ggs <- purrr::map(doublet_vars, function(x) {
p <- DimPlot(tonsil, group.by = x, pt.size = 0.3) +
ggtitle(x) +
theme(
plot.title = element_text(size = 14, hjust = 0.5),
axis.title = element_text(size = 13),
axis.text = element_text(size = 11),
legend.text = element_text(size = 11)
)
p
})
names(doublet_ggs) <- doublet_vars
doublet_ggs
## $HTO_classification.global
##
## $scrublet_predicted_doublet
##
## $has_high_lib_size
# Zoom in for cell hashing
umap_hashing <- tonsil@reductions$umap@cell.embeddings %>%
as.data.frame() %>%
mutate(HTO_classification.global = tonsil$HTO_classification.global) %>%
ggplot(aes(UMAP_1, UMAP_2, color = HTO_classification.global)) +
geom_point(size = 0.01, alpha = 0.5) +
facet_wrap(~HTO_classification.global) +
theme_classic()
umap_hashing
canonical_markers <- c("CD79A", "CD79B", "CD3D", "CD3E", "NKG7", "LYZ", "FDCSP")
canonical_markers_gg <- purrr::map(canonical_markers, function(x) {
p <- feature_plot_doublets(seurat_obj = tonsil, feature = x)
p
})
names(canonical_markers_gg) <- canonical_markers
canonical_markers_gg
## $CD79A
##
## $CD79B
##
## $CD3D
##
## $CD3E
##
## $NKG7
##
## $LYZ
##
## $FDCSP
We want to be particularly careful with two cell types: pre-plasmablasts and tingible-body macrophages, since they might be confused by doublets (specially using scrublet). Thus, we will visualize a prePB and apoptosis signature:
# Define signatures
apoptosis_signature <- AnnotationDbi::select(
org.Hs.eg.db,
keys = "GO:0006915",
keytype = "GO",
columns = "SYMBOL"
)$SYMBOL
apoptosis_signature <- unique(apoptosis_signature)
prePB_signature <- c("RASSF6", "FRZB", "HOPX", "BRNL9", "FGFR1")
signatures_list <- list(
prePB_signature = prePB_signature,
apoptosis_signature = apoptosis_signature
)
# Compute score
tonsil <- AddModuleScore(
tonsil,
features = signatures_list,
name = c("prePB_signature", "apoptosis_signature")
)
# Visualize
prePB_gg <- feature_plot_doublets(
seurat_obj = tonsil,
feature = "prePB_signature1"
)
apopotosis_gg <- feature_plot_doublets(
seurat_obj = tonsil,
feature = "apoptosis_signature2"
)
prePB_gg
apopotosis_gg
tonsil
## An object of class Seurat
## 28504 features across 347262 samples within 1 assay
## Active assay: RNA (28504 features, 3000 variable features)
## 3 dimensional reductions calculated: pca, umap, harmony
table(tonsil$HTO_classification.global)
##
## Doublet NA Singlet
## 80577 50072 216613
tonsil <- subset(tonsil, subset = HTO_classification.global != "Doublet")
tonsil
## An object of class Seurat
## 28504 features across 266685 samples within 1 assay
## Active assay: RNA (28504 features, 3000 variable features)
## 3 dimensional reductions calculated: pca, umap, harmony
saveRDS(tonsil, path_to_save)
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
##
## Matrix products: default
## BLAS: /apps/R/3.6.0/lib64/R/lib/libRblas.so
## LAPACK: /home/devel/rmassoni/anaconda3/lib/libmkl_rt.so
##
## locale:
## [1] LC_CTYPE=C LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.4 purrr_0.3.4 readr_1.3.1 tidyr_1.1.0 tibble_3.0.1 ggplot2_3.3.0 tidyverse_1.3.0 Matrix_1.2-18 Seurat_3.2.0 org.Hs.eg.db_3.10.0 AnnotationDbi_1.48.0 IRanges_2.19.14 S4Vectors_0.24.4 Biobase_2.45.1 BiocGenerics_0.32.0 BiocStyle_2.14.4
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.7 plyr_1.8.6 igraph_1.2.5 lazyeval_0.2.2 splines_3.6.0 listenv_0.8.0 digest_0.6.20 htmltools_0.4.0 viridis_0.5.1 fansi_0.4.1 magrittr_1.5 memoise_1.1.0 tensor_1.5 cluster_2.1.0 ROCR_1.0-11 globals_0.12.5 modelr_0.1.8 colorspace_1.4-1 blob_1.2.1 rvest_0.3.5 rappdirs_0.3.1 ggrepel_0.8.2 haven_2.3.1 xfun_0.14 crayon_1.3.4 jsonlite_1.7.2 spatstat_1.64-1 spatstat.data_1.4-3 survival_3.1-12 zoo_1.8-8 ape_5.3 glue_1.4.1 polyclip_1.10-0 gtable_0.3.0 leiden_0.3.3 future.apply_1.5.0 abind_1.4-5 scales_1.1.1 DBI_1.1.0 miniUI_0.1.1.1 Rcpp_1.0.6 viridisLite_0.3.0 xtable_1.8-4 reticulate_1.16 bit_1.1-15.2 rsvd_1.0.3 htmlwidgets_1.5.1 httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.1 ica_1.0-2 farver_2.0.3 pkgconfig_2.0.3
## [55] uwot_0.1.8 dbplyr_1.4.4 deldir_0.1-25 here_0.1 labeling_0.3 tidyselect_1.1.0 rlang_0.4.10 reshape2_1.4.4 later_1.0.0 munsell_0.5.0 cellranger_1.1.0 tools_3.6.0 cli_2.0.2 generics_0.0.2 RSQLite_2.1.2 broom_0.5.6 ggridges_0.5.2 evaluate_0.14 fastmap_1.0.1 yaml_2.2.1 goftest_1.2-2 knitr_1.28 bit64_0.9-7 fs_1.4.1 fitdistrplus_1.1-1 RANN_2.6.1 pbapply_1.4-2 future_1.17.0 nlme_3.1-148 mime_0.9 xml2_1.3.2 compiler_3.6.0 rstudioapi_0.11 plotly_4.9.2.1 png_0.1-7 spatstat.utils_1.17-0 reprex_0.3.0 stringi_1.4.6 lattice_0.20-41 vctrs_0.3.6 pillar_1.4.4 lifecycle_0.2.0 BiocManager_1.30.10 lmtest_0.9-37 RcppAnnoy_0.0.16 data.table_1.12.8 cowplot_1.0.0 irlba_2.3.3 httpuv_1.5.3.1 patchwork_1.0.0 R6_2.4.1 bookdown_0.19 promises_1.1.0 KernSmooth_2.23-17
## [109] gridExtra_2.3 codetools_0.2-16 MASS_7.3-51.6 assertthat_0.2.1 rprojroot_1.3-2 withr_2.4.1 sctransform_0.2.1 mgcv_1.8-31 hms_0.5.3 grid_3.6.0 rpart_4.1-15 rmarkdown_2.2 Rtsne_0.15 shiny_1.4.0.2 lubridate_1.7.8